General recursive solution for one-dimensional quantum potentials: a simple tool for 
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A revision of the recursive method proposed by S.A. Shakir [Am. J. Phys. 52, 845 (1984)] to solve 
bound eigenvalues of the Schrodinger equation is presented. Equations are further simplified and 
generalized for computing wave functions of any given one-dimensional potential, providing accurate 
solutions not only for bound states but also for scattering and resonant states, as demonstrated here 
for a few examples. 
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I. INTRODUCTION 

Predicting physical properties of quantum systems has 
basically been treated as an eigenvalue problem. Since 
from single atoms and molecules up to nanostructured 
devices, chemical and electronic properties are deter- 
mined by finding the eigenstates of the system's Hamilto- 
nian. Methods of solving eigenvalue/eigenfunction prob- 
lems are therefore an inevitable part of the quantum 
mechanic theory. In real systems, even one dimen- 
sional problems can be quite difficult to solve analyti- 
cally and often approximative, such as in perturbation 
theory, and/or numerical methods are requiredi^^^ 
Semiconductor heterostructures are good examples of one 
dimensional systems of technological interest where nu- 
merical solutions of the Schrodinger equation (SE) are 
often employed for self-consistent studies of electronic 
devices.— i 8 ' 9 ' 10 ! 11 Studies of vibrational bound states of 
molecules in physical chemistry, as well as of other 
oscillatory systems in atomic and nuclear physics are 
treated as one dimensional problems when summarized 
in the radial SE; they consist another class of impor- 
tant examples in which numerical methods are also often 
employed, 12 ' 13 ! 14 . 15 

Although many procedures and algorithms are avail- 
able in nowadays for solving the SE numerically — see for 
instance Refs. [HI, Qj], and references therein — they are 
suitable to address specific problems on distinct fields of 
applied and theoretical physics. There is still a need for a 
general treatment that can be summarized into a single 
routine capable to solve equally well problems of con- 
fined modes by arbitrary potential as well as scattering 
problems in open systems. Moreover, all available numer- 
ical methods are essentially mathematical solutions of a 
differential equation, and hence the level of expertise re- 
quired for dealing with these methods has kept them out 
of most physics and chemistry classrooms, even on grad- 
uate courses, where the demonstrative examples are still 
limited to a few cases of quantum potentials for which 
the analytical solutions are known. 

More than 20 years ago, S.A. Shakir— presented a re- 
cursive method to find bound eigenvalues of the SE. It 
was a relatively simple method, closely related to those 
used in the field of optics for calculating the Fresnel am- 
plitudes. In despite of its simplicity, and even of the fact 
of have been rediscovered a few years latter ; 19 ' 20 ! 21 this 



method has passed unnoted from the applied research 
fields until very recently^ 

In this work, instead of concerning with mathemati- 
cal methods to solve the SE for bound and scattering 
states, we are concerned with alternative methods capa- 
ble of providing better physical insight about the quan- 
tum world. In other words, methods that provide solu- 
tions of the SE without having to solve it mathematically. 
One class of alternative methods, which we call physical 
methods, is possible in one-dimension. They are based 
on the elemental action of a given stationary force field 
on the particle's wave function. Action that occurs at ev- 
ery step of a discretized potential energy function and it 
is qualitatively always the same, reflecting and transmit- 
ting the incident wave as in the Shakir's method. 18 How- 
ever, here, the original formalism is simplified and gener- 
alized into a procedure for computing wave functions of 
any given potential whatever related to scattering, reso- 
nance, tunneling, or bound states, as demonstrated for 
a few examples. Physical method at their actual stage 
of development might require higher computational cost 
than mathematical methods, but offers versatility and 
exactness in return, as well as the fact that they can be 
explained and used even by undergraduate students and 
no-experts in numerical methods. 



II. THEORY 

For a given potential energy function U (x), it is always 
possible to define the interval X = [xq, xn] such that out- 
side its boundaries, i.e. for x X, either the potential 
is constant or the amplitude of the wave function van- 
ish before appreciable variation of the potential. In this 
interval, the potential has been discretized into N steps, 
and we seek for stationary states with wave functions in 
the form 



3=0 



in which 



(1) 



(2) 



and 5j(x) = 1 for 1 and otherwise. The 

under arrow indicates that this format will be used here 
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FIG. 1: Discretization of a energy potential function U(x) 
(dashed line), providing reflection, Rj and Rj, and transmis- 
sion, Tj and Tj, coefficients at every step of the discretized 
function (solid line). Coefficients with and without bar stand 
for right-hand and left-hand scattering solutions, respectively. 
Wavevectors kj at each interval of constant potential value 
U(xj) is also shown. 

for left-hand solutions, as better explained below. 

Incremental distances dxj — Xj-\~\ — Xj , 3X6 small 
enough to guarantee a good numerical solution when tak- 
ing U{x) = U(xj)Sj(x), as shown in Fig. 1, so that 

kj = tp^J E - U( Xj ) (3) 



where if — y 2m/ h 2 . E and m stand for energy and mass 
of the particle, respectively. 

To calculate the Aj and Bj amplitudes recursively, ini- 
tial values at the boundaries of the interval X have to 
be provided. For instance, normalized solutions for inci- 
dence from the left-hand side, i.e. incidence from x < xq, 
must have Bn = and Aq = 1. Conservation of the 
probability current imposes that 

i^j-i = and i^j-i = tfj ( 4 ) 

at every step of the potential, i.e. at every Xj. The 
prime symbol indicates the first derivative in x. Left- 
hand solutions are obtained by applying these conditions, 
Eqs. (J4j> , first at xn where Bn = and 

A N = An-! . 2kN ;\ e ik »-^ x "- x »-^ = A N -iT N , 
k N -i + k N 

(5) 

and then at xn-i where 

R A kff-1 - KN 2 ik N -l(x N -X JV-l) — A R 

Bn-i—An-i- — — e — An-iHn 

kn-i + k N 

(6) 

and 



A a 2fcjy_ 2 ifcjv_2(a;jv-i— xn-z) A T tT\ 

An-x-An-2-tt j Vr — 77u TT \ e '-An-2-Ln-i, (7) 

(fcjV-2 - kN-l)riN + {kN-2 + KjV-lJ 



and then at xn-2 where 



R _ A (fcjV-2 + k N -l)RN + (fcjV-2 ~ fcjV-l) ^k N - 2 (x N - 1 -x N - 2 ) A R fo\ 

Bn-2—An-2J7 7 \TJ —77 TT " e V ' - A N -2Rn-1 (=>) 

{kN-2 - k N -i)riN + {kN-2 + Kn-1) 



r 



and An-2 is analogous to Eq. ([7]) when replacing N by 
7V-1. 

This is a recursive procedure, repeatable since j = N 
down to j — 1 as summarized by 

Aj = A, ,7, and /,', = A j R j+1 (9) 

where 

3 {kj-! - kj)R j+1 + (kj-x + kj) 

(10a) 

and 

R . = fe-i + kj)R j+1 + (fcj-i - k 3 ) c2ik ^ i(x ^_ Xi i) 
3 (kj-! - kj)R j+1 + {kj^x + kj) 

(10b) 

All amplitudes are proportional to Aq, and within 
the condition assumed when defining the interval X, 
Rn+i = 0. 



Solutions for incident particles from the right-hand 
side, i.e. incidence from x > xn, are obtained in analo- 
gous procedure, but for the sake of simplification in the 
final format of the recursive equations it is convenient to 
redefine the amplitudes as follow 

Aj = C j+1 e- ik i dx * and //, D j+1 t •"' dx . (11) 

and hence the wave function ^(x) of the stationary states 
will be calculated by using 

ipj(x) = C j+1 e ik ^ x - x ^ + D j+1 e- ik ^ x - x ^ (12) 

instead of ipj(x). 

By applying the continuity of ip j and ip'j, as in 

Eqs. (HI), analogous deduction of that in Eqs. (©-© can 
be carried out but now starting from xi where C\ = 
up to xn+i where Dn+i = 1 (for normalized solutions). 
It leads to the recursive equations 

Bj = D j+1 fj and C 3 = DjRj-x (13) 
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of the right-hand solutions where 



T, 



— ' ■ 



and 



(14a) 



(kj + kj-tjRj-i + (kj - fcj-i) 2ifc ,(x, +1 -x,) 



(kj — k^-ijRj-i + (kj + kj- 







(14b) 

for j = 1,2, ...,JV and R = 0. 

At certain energies E = £ v , quantum confinements 
can occur at local minima of the potential in which 
U(x) < £„; given that each minimum is bounded on 
both sides by non tunnelable barriers. Around one of 
such local minima, stationary wave solutions only ex- 
ist when the waves reflected at both bounding "walls" 
interfere constructively to each other at any instant of 
time. It means that, left-hand and right-hand scatter- 
ing solutions must coexist in the confinement region, i.e. 
ip j(x) = ip j(x), which takes us back to Eqs. (fTTj) . Since 

Bj = AjRj + i and Cj + i = Dj +i Rj, as given by Eqs. © 
and (fl~3|) . all amplitudes cancel each other out in Eq. (fTTj) 
so that RjRj+i = exp(2ikjdxj). This latter relationship 
provides a criterion for numerical determination of the 
allowed modes in any confinement region of the potential 
since a minimization function such as 



f(E) 



RjR 



(15) 
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FIG. 2: Resonant tunneling in a double barrier with V-shaped 
well at the middle, (a) Transmission coefficient T, Eq. (JT7J, 
as a function of the particle's energy E; plot resolution is 
dE = 0.002eV. (b) Potential energy function of the barrier 
and wave functions for the two energies (0.406eV and 0.45eV) 
pointed out by arrows in (a). A set of curves is shown for 
each wave function solution, corresponding to ^(x)e _I ^' for 
several values of ujt. The amplitude of the wave function 
for E = 0.45eV has been magnified by a factor of 2 regard- 
ing the amplitude of the resonant wave. Vertical displace- 
ments of both sets of curves are not related to the energy 
scale at left. Solutions computed for m = 511keV/c 2 , i.e. 
ip = 5.1232eV -1 / 2 nm -1 in Eq. and an uniform step width 
of dxj — 0.02nm in the interval a;jv,o = ±10nm (N=1000). 



can be defined with j running over all values where 
U(xj) < E. Then, the condition f(E) = is fulfilled 
only for a set of discrete £ v values of energy, correspond- 
ing to the eigenvalues of the one dimensional SE at a 
given local minimum of U(x). 

Eigenfunctions of the energy eigenvalues £ v , are 
computed as either ip j( x ) 01 j( x ) solutions in the clas- 
sical allowed region around the minimum. But both so- 
lutions are required to compute the evanescent parts of 
the eigenfunctions inside the bounding walls, and then it 
is necessary to match the amplitudes of these solutions 
at some point. If the hth step, for which U(xf l ) < £ Vl 

is taken as the matching point, -0 ft = ^ h and hence 

— > < — 

D h = A h (l + R h+1 )/(l + Rh-i). By providing an arbi- 
trary value to Ah, such as Ah = 1, the non-normalized 
eigenfunction 

f i/jj(x), if h<j< N 
*„(*) = !>(*) x | £ j(jB)) i£1 < j<h (16) 

can be calculated over the entire interval X. Since 
S = X^jlJo 1 \^v(xj)\ 2 Axj has to be equal to unit for a 
normalized eigenfunction, y/S is the normalization fac- 
tor. 



III. EXAMPLES AND DISCUSSIONS 

A simple application of this recursive method is in find- 
ing transmission coefficients, 

T=\A N /A \ 2 = 1-\B Q /A \ 2 , (17) 

across arbitrary potential barriers. It is known as be- 
ing able to provide exact results even for truncated po- 
tentials where the most common approximative methods 
are used to fail, as compared elsewhere^ Besides ex- 
actness, another advantage of the present formalism is 
that it is not just computing the values of T, but the 
entire wave function instead. Hence, for more complex 
potential barriers, as those usually find in heterostruc- 
tures and nanodevices where resonance can take place 
inside the barrieS ) 17 i 22 ' 24 i 25 amplitudes and evanescence 
times of resonant modes are simultaneously computed. 

As an example, consider the potential given in Fig. 2, 
which is a double barrier with V-shaped well the middle. 
By scanning the energy values E, Eq. ([3]), in a chosen 
range from E m i n to E max in Ne steps of width dE so that 
E n = E min + (n - l)dE and E Ne = E max , a T(E) curve 
is obtained as the one shown in Fig. 2(a). By storing 
the coefficients Aj(E n ) and Bj(E n ) when calculating the 
T values, high-resolution plots of the wave functions in 
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the interval X also as a function of energy are already 
available, i.e. 

N 

*(x, E n ) ~ J2 S j( x )iM E n) + B 3 (E n )}, (18) 

3=0 

given that dxj << Xj = 27r/|fcj|. Hence, amplitudes 
and attenuations of the waves through the entire barrier 
can be visualized. For this particular potential, it shows 
that the mode undergoing resonant-tunneling, at E — 
0.406eV, has two maxima in the well as can be seen in 
Fig. 2(b). 

Optionally one can compute only the Rj's and X^'s, 
Eq. (US}, to obtain T = UI^i^l 2 and either stop the 
calculation there if only the T(E) curve is desired or com- 
pute the Aj's and Bj's, Eqs. ([9]), later for a chosen set of 
energies. The options on how the present formalism can 
be programmed for optimizing computational tasks are 
detailed in Appendix A. 

Although the stored set of solutions, Eq. (fl~8|) , could 
be used to study the scattering of wave packets in the 
barriers, obtaining normalized wave packets would not be 
as straightforward as possible because the solutions are 
evenly spaced in energy, not in momentum. Therefore, 
if one wants to compute wave solutions only once, and 
also study wave packets composed of these solutions as 
well as their transmission coefficients, it is better to first 
define the form of the initial wave packet. 

Gaussian wave packets has become a benchmark in 
computing time evolution of this sortfii hence let use 

$( Xj 0) = 1 P i^ p -(x-x a fl2al 

as reference for our initial free-particle wave packet in 
t = 0, centered at x , and described by a single mode 
of energy Eq so that kq — (p\/E~Q. It is normalized 
since J \^(x, 0)\ 2 dx = 1, and its full width at half 
maximum (FWHM), W, can be chosen by providing 
(j x = W7V2hi4 ~ 0.6W. 

Numcrov-type of solutions of the time-dependent SE 
can handled the time evolution of such wave packet, 
Eq. (TiT)|) , through a given potential barrier jil but in the 
present recursive method it is not possible to start with 
such expression of ^(x,0) since it is an artificial expres- 
sion that happens to have the same shape, in t = 0, of 
the actual wave packet 

, We 

ttfo t) = -±= V c n *(x, E n )e- lE - l / n (20) 

where 

C n = 1 e -(Kn-K0) 2 /2al (21) 

and cr/j = a~ x . For sake of normalization, the wave vec- 
tors K n must be equally spaced, i.e. K n +i — K n = dn, so 




FIG. 3: Scattering of a gaussian wave packet in a double bar- 
rier with V-shaped well (thick line, Fig. 2). Time evolution 
of ^{x,t) (thin line) and |$(a;,t)| 2 (gray shaded curves), in- 
creasing from (a) to (d), as indicated aside of each snapshot. 
Wave packet obtained according to Eqs. (|20p and (|21l) . as de- 
scribed in the text. Group velocity vq = 0.378nm/fs (kq = 
3.264nm _1 ), initial position at xo = — 60nm from the center 
of the barriers, and FWHM of 25nm (a x = 1/<t& = 15nm). 

that 

N E 

Y J \c n \ 2 dn = l. (22) 

n=l 

Since the c n 's out of the range K n — kq — ±3.5crfc are of 
negligible amplitude, the Ne modes with wave vectors 
K\, K2, ■ ■ ■ , kn e are taken in this range to compose the 
wave packet, whose spectrum of energy is given by E n = 
(n n /ip) 2 . Energy difference between adjacent modes, 

dE n = E n+ i - E n ~ 2n n dK/ip 2 , 

increases as n runs toward Ne- Hence, the maximum 
interval of time by which the evolution of the wave packet 
can be studied is 

h s h ( AE\ 
At < — w (Ne - 1)— ;— 1 

" 2(E Ne -E Ne _ 1 ) v 4AE \ 2E J 

(23) 

when the energy range Eq ± AE comprises all E n 's. 

To be more illustrative, let construct an initial guassian 
wave packet within a previously chosen range of energy, 
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FIG. 4: (a) Vibrational energy eigenvalues estimated via 
minimization of f(E), Eq. (|15|l . for a Lennard- Jones type 
of potential and for two rotational states, with J = 
(solid line) and with J = 8 (dash-dot line). Step width 
dxj = 0.95pm. Potential energy functions given by 
U(x) = A/x 12 - B/x 6 + J(J + 1) /{ipxf where A = 0.124 x 
10" 12 eV ■ nm 12 , B = 1.488 x 10~ 6 eV • nm 6 and m = 
469.4MeV/c 2 . (b) Eigenfunctions of the three first vibrational 
states for the case J = 0, indicated by arrows in (a), are shown 
in the same scheme as before (Fig. 2), i.e. in sets of 20 curves 
per time period of amplitude oscillation. Wave amplitudes 
are in a common scale. Plot of U(x) (thick-solid line), J = 0, 
is also shown. 



for instance, with E n in the narrow range Eq ± 0.058eV 
around the energy Eq = 0.406eV of the resonant mode 
in Fig. 2. It leads to a k ~ ipAE/7y r E ~ = O^enm" 1 
and to W = 25.0nm as the FWHM of |f(x,0)| 2 . By 
composing the wave packet with 101 modes, our limit of 
time for studying its evolution is At < 1654fs, Eq. (|2"3")) . 
Fig. 3 shows the scattering process of such wave packet 
by the double barrier (Fig. 2) as computed via Eq. ([20]) . 
The reflection of the entire wave packet occurs in a time 
interval no longer than 170fs, Fig. 3(a) through Fig. 3(c), 
but the resonant mode trapped in the well survives for 
much longer than that, Fig. 3(d). Its evanescence in time 
is dictated by an exponential decay of time constant of 
184. 5fs regarding the probability of finding the particle 
inside the well after t = 232. 3fs, Fig. 3(c). 

In case of potential with localized minima not reach- 
able by tunneling from either sides, computing both left- 



hand and right-hand solutions is necessary. The proce- 
dure is very similar to the previous one used to obtain the 
ty(x,En) set of wave functions evenly spaced in energy 
where E n+ i — E n = dE. But, most of these unilateral 
solutions do not exist under confinement. Hence, prior 
to calculate Ay's, Bj's, and Tj's from Eqs. §§§ and (jlOap . 
we first compute the Rj's and -Rj's, Eqs. (|10b[) and (|14b[) . 
to generate the f(E) plot according to Eq. (fTSj) . In Fig. 
4(a) there are examples of such plot for a Lennard- Jones 
potential, a type of potential used to describe diatomic 
molecules and that has a single minimum. Only after 
providing an energy eigenvalue, £ v , all other required co- 
efficients for plotting the eigenfunction ^ v (x), Eq. (fl6|) . 
are calculated. For u = 1,2, and 3, the eigenfunctions are 
shown in Fig. 4(b) where the one-dimensional variable x 
stands for the nuclei distance in the molecule. Although 
reduced mass, equilibrium distance, and ionization en- 
ergy values are close of those for the H2 molecule, the 
width of the depression in U (x) is much narrow than in 
the actual molecule^ It means a stronger restoring force 
and hence larger energy gaps in between adjacent vibra- 
tional states, which is used here for illustrative purposes 
only. In this hypothetical molecule, the non-constancy 
of the gap values is more evident, as can be seen in Fig. 
4(a). It indicates how the exact solution would differ 
from the harmonic oscillator approximation. 

Computational task involved in finding the energy 
eigenvalues E = £ v in which f(E) = can be enor- 
mous depending on number and accuracy of the desired 
eigenvalues. There will be no difficult in determining 
a few eigenvalues with reasonable accuracy, for instance 
£ 1 = -3.4416eV, £ 2 = -1.8793eV, £ 3 = -0.8660eV, 
and £4 = — 0.2991eV in Fig. 4(a), are obtained at once 
in just a few seconds and with an accuracy of dE/2 = 
0.0022eV (N E = 1000); x = 0.002nm, x N = 0.2nm, 
and dxj = 0.0005nm (N = 396). The problem relies in 
finding many eigenvalues, with high accuracy, occurring 
at different densities along a broad energy range. But, 
practical needs of finding all possible eigenvalues are un- 
usual, although it can be realized by developing search 
routines to reduce computational time since the method 
does not have any intrinsic limitation; as for instance the 
frequency distortion, or phase-lag, errors of the Numerov- 
type methodSi 12 i 13 

Potential functions with singularities do not compro- 
mise the applicability of the presented method since only 
solutions with null wave functions at the singularities are 
those with physical meaning. In other words, when rep- 
resenting the potential by a discretized function, a finite 
value has to be assigned to the singularity. It compro- 
mises mostly the modes with non-null wave functions at 
the singularity, but these are artificial modes that do not 
exist in the actual physical system. Occurrence of artifi- 
cial modes has been observed when solving the potential 
JJ(x) = -e 2 /[M + |e|] with e -> 0. In this case, the 
odd solutions with respect to the singularity provide the 
well-known eigenvalues of the hydrogen atom. 

Another situation that can also be addressed by the 
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FIG. 5: Energies and stationary wave functions of con- 
fined modes in a potential with two minima. (a) f(E) 
curves for the right (black line) and left (gray line) min- 
ima, showing the allowed energies in the wells. (b), 
(c), (d) Confined modes at the energies — 62.68(12)meV, 
-57.27(13)meV, and -51.86(12)meV, as indicated by ar- 
rows in (a). 9(x)e~'' u ' t (thin lines) for several values of 
tut, and ^(x)] 2 (gray shaded curves) are shown as well as 
the plot of U(x) = A(x - a) 2 + Bexp{-(x - a - 6) 2 /a 2 } + C 
(thick line); A = 4.0meV (= 2.4meV) for x < a (x > a), 
B = 450meV, C = -500meV, 8 = 0.5nm, and a = lOnm. 
xn,o = a, ± 20nm and dxj = O.lnm (N=400). Particle's mass 
m = 1.022MeV/c 2 . 



recursive method occurs when the potential function has 
closely spaced minima, as for instance in the potential 
shown in Fig. 5. In the present formalism, the two min- 
ima (or wells) of this potential can be treated separately 
after establishing distinct intervals in X for calculating 
the function f(E). For each E n g [E min , E max ], we first 
compute the Rj's and Rj's in the entire interval X that 
contains the wells, i.e. from xo to xn- Then, by se- 
lecting a position xt somewhere in between the wells as 
indicated in Fig. 5(b), two f{E n ) values are obtained: 
one for Xj <E [xo,Xb] and another for Xj € [xb+i, £;v], and 
always obeying that U(xj) < E n . Both f(E) curves thus 
obtained are shown in Fig. 5(a). As tunneling across the 
barrier separating the wells become relevant, the minima 
of the f(E) curves are coincidental. It means that the 
energy levels, or eigenvalues £„, in each well are indepen- 
dent from each other only if tunneling is negligible. It 
is more evident when visualizing the wave functions of 
the confined modes shown in Figs. 5(b), 5(c), and (d). In 
practice, potential with two minima are found for Cooper 
pairs (2 spin-coupled electrons) in josephson-junction of 
superconducting materials — 

In semiconductor heterostructures, the validity of this 
recursive method in its present formalism holds true 
when the charge carrier's effective mass can be ap- 
proximated by a constant value through out the struc- 
ture. However, accounting for effective mass variation 
across heterojunctions has only been possible in cases of 



abrupt interfaces and by using appropriated boundary 
conditions! 28 ' 29 

Computing wave functions recursively is limited to 
problems treatable in one dimension, such as spherically 
symmetric potentials, or cases where the potentials in 
orthogonal directions are independent from each other, 
as for instance U(x,y) = Ux(x) + Uy(v), in which the 
wave functions are obtained by the product of indepen- 
dent one-dimensional solutions. In all other cases of two 
and three dimensions, it is not possible to extend this 
recursive method since the number of unknown ampli- 
tude coefficients are more than the number of equations 
written by matching the wave function and its derivative 
at every steps of the discretized potentials. Hence only 
mathematical methods for solving 2D and 3D differential 
equations are available ! 3 ' 5 ' 29 ' 30 



IV. CONCLUSIONS 

This work has demonstrated that any one-dimensional 
quantum potential, wherever representing an open or 
closed system, can be solved by a single method and 
within good numerical accuracy. It is a simple tool ready 
to be used in applied physics or just in developing illus- 
trative examples of quantum mechanics. It should also 
be emphasized that, in the presented formulation, the 
quantum potentials are solved by a physical treatment 
where the particle's wave function is molded by the same 
elemental action of the given force field. This is a com- 
pletely different approach of the one that is usually taken 
where, once the potential is known, the problem became 
the mathematical solution of a differential equation. 



APPENDIX A: ROUTINES FOR COMPUTING 
WAVE SOLUTIONS RECURSIVELY 

There are several ways to program the recursive for- 
malism presented in this work, depending on what infor- 
mation is desired. The flowchart in Fig. 6 describes the 
basic structure of the routines used to elaborate the ex- 
amples in Figs. 2 to 5. Once the U{x) function and the 
discretization of X is defined, the user can choose either 
route © or @ to scan a given energy range in fixed incre- 
ments of energy, dE, or wavevector, dn, respectively. The 
latter route is suitable to study time evolution of wave 
packets by ending through route ©, i.e. ©— >®^©^@. 
Eigenvalues {£\, £2, ■ ■ •}, of confined modes are obtained 
from the f(E) curve in route @, and used as input to 
provide the eigenfunctions ^ u (x) through route ®. The 
straight route @— >© are exclusive for users interested 
only in transmission coefficients across arbitrary poten- 
tial barriers. 
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FIG. 6: Flowchart of possible routines for studying one dimensional potentials with the present recursive formalism. Trans- 
mission coefficients, T(E) curves only: @^©^©. Wave function ^(x,E) plots, and/or T(E) curves: ©^©^©. Time 
evolution of wave packet ty(x,t) plots: @^©^©^@. Eigenvalues S v , and eigenfunction ^ v (x) plots: ®— >@— >®. 
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